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Abstract 

In the literature, two quite different phase-field formulations for the problem of alloy solidification 
can be found. In the first, the material in the diffuse interfaces is assumed to be in an intermediate 
state between solid and liquid, with a unique local composition. In the second, the interface is 
seen as a mixture of two phases that each retain their macroscopic properties, and a separate 
concentration field for each phase is introduced. It is shown here that both types of models can be 
obtained by the standard variational procedure if a grand-potential functional is used as a starting 
point instead of a free-energy functional. The dynamical variable is then the chemical potential 
instead of the composition. In this framework, a complete analogy with phase-field models for the 
solidification of a pure substance can be established. This analogy is then exploited to formulate 
quantitative phase- field models for alloys with arbitrary phase diagrams. The precision of the 
method is illustrated by numerical simulations with varying interface thickness. 

PACS numbers: 64.70.Dv,81.30.Fb,05.70.Ln,81.10.Aj 
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I. INTRODUCTION 



The development of the phase-field method has led to tremendous progress in the mod- 
elling of pattern formation during solidification, due to its capability to simulate complex 
;ime-dependent and three-dimensional morphologies with relatively simple numerical codes 
l|-|4(. The general principle of this method is to describe a multi-phase system by a set 
of phase fields which take constant values in each of the bulk phases and vary smoothly 
through interfaces of a characteristic thickness W. The equations of motion for the phase 
fields and their coupling to the local thermodynamic state variables (temperature, density, 
composition etc.) can be obtained, following the basic principles of out-of-equilibrium ther- 
modynamics, by taking a variational derivative of a free energy functional, which is generally 
of the Ginzburg-Landau type. Mean-field approximations can be used to relate the param- 
eters that appear in this functional to microscopic quantities, and the phase fields can often 
be interpreted as order parameters. 

Generally, the equations that result from the straightforward application of these prin- 
ciples are not suitable for obtaining quantitatively accurate simulation results on solidifica- 
tion microstructures. The reason is that the characteristic natural thickness of the diffuse 
solid-liquid interfaces is a few times the interatomic distance, whereas solidification patterns 
typically exhibit length scales ranging from 1 to 100 /xm. Even with the help of modern 
computers and multi-scale algorithms, both of these scales cannot be resolved at the same 
time. Therefore, in order to simulate solidification microstructures, the thickness of the dif- 
fuse interfaces in the phase-field model has to be artificially enlarged, sometimes by two or 
three orders of magnitude. Quantitative results can only be expected if both the equilibrium 
and kinetic properties of the interfaces remain unchanged under this procedure. 

To achieve this goal, it is helpful to adopt a phenomenological point of view: the phase 
field is seen as a smoothed indicator function (as opposed to a physical order parameter or 
density), and all equilibrium quantities and transport coefficients are interpolated between 
the phases with smooth functions of the phase fields that can be freely chosen. This freedom 
can be exploited to construct phase-field models with special properties. In particular, 
rescaling of the interface thickness is greatly simplified in models where bulk thermodynamics 
and interfacial properties can be controlled separately. 

The phase-field models for the solidification of a pure substance published in the literature 



are all quite similar [l|, |5K7j , in the sense that they use the same set of fundamental fields (a 
phase field and the temperature field), and that the structure of the equations is the same. 
One reason for this universality is that a simple and intuitive formulation of the model in 
terms of these fields yields indeed, as will be detailed below, a model in which this separation 
is achieved. Therefore, the development of the "thin interface limit" l|, |8[ which has paved 
the way to quantitative simulations of dendrites l|, |9[ did not require a change in the model 
formulation. 

The situation is more complicated for alloy phase-field models. Two different approaches 
with quite distinct philosophies have been pursued in parallel. The first, which will be called 
the "coarse-graining" approach in the following, generalizes the pure substance model by 
introducing a concentration field in addition to the temperature field and by writing down 
a free e nergy 



tration 
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'unctional that depends on the phase field, the temperature, and the concen- 
121 ]. The local values of these three fields are - in principle - the coarse-grained 
counterparts of the microscopic structural order parameter, temperature, and concentration 
fields, and the interface is seen as a narrow region in space where all of these quantities 
can exhibit rapid spatial variations. In contrast, the second approach, called "two-phase" 
approach in the following, treats the interfaces as a mixture of two phases, each of which 
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141 ] . In this approach, the phase 



retains its bulk properties even inside the interface 
field represents the volume fraction of one of the phases. Moreover, a separate concentration 
field is defined for each phase, and the physical concentration field is obtained by a weighted 
average depending on the value of the phase field. The introduction of two separate con- 
centration fields adds a supplementary degree of freedom which has to be removed from the 
problem; this is done either by a specific partition relation [l3| or by the condition of local 
equilibrium between the coexisting phases |l4j . 

In the models of the "coarse-graining" type, the model structure generally leads to an 
intrinsic coupling between bulk and interface properties, which makes simulation results 
dependent on the chosen interface thickness. Only for the case of dilute binary alloys, a 
specific interpolation of the thermodynamic properties through the interface has been devel- 
oped 15|, Il6] which overcomes this constraint and makes quantitative simulations possible. 
In the two-phase approach, which is more phenomono logical from the outset, bulk and inter- 
facial properties are decoupled by construction. However, the removal of the extra degree of 
freedom introduced by the model formulation generally requires the solution of a nonlinear 



equation in each point of the interface, and thus adds significant computational complexity. 

The purpose of the present paper is to show that the coarse-graining approach can be 
easily extended to more complex alloy systems if instead of a free energy functional a grand- 
potential functional is used to generate the equations of motion. Furthermore, an analysis of 
the resulting model shows that it is in fact perfectly equivalent to the two-phase model, which 
offers the possibility to reinterpret and simplify the latter. The fundamentals underlying 
these findings can be stated quite simply. The motion of interfaces is controlled by the 
transport of a conserved extensive quantity: energy for a pure substance, and chemical 
species for isothermal alloy solidification. In sharp-interface models, this fact is expressed by 
two separate laws: a transport equation in the volume and a conservation law of Stefan type 
at the moving boundary. In contrast, the two-phase equilibrium at interfaces is controlled 
by the intensive quantity that is conjugate to the conserved one: temperature for pure 
substances and chemical potential for alloys. It turns out that the pure substance model 
has been formulated from the start in terms of a phase field and the intensive variable 
(temperature), whereas alloy phase field models are traditionally formulated in terms of a 
phase field and the composition, which is a density of the extensive variable (number of 
solute atoms). To obtain a model for alloys that has properties analogous to those of the 
pure substance model, it is sufficient to choose a formulation in terms of a phase field and 
the chemical potential, and to switch to the appropriate thermodynamic potential, which is 
the grand potential. Of course, the fact that the quantity analogous to the temperature is 



the chemical potential is well known 



171 ] and has been extensively used in sharp-interface 



models as well as in a few specific phase-field models [18l-l20|. but its general consequences 
have so far not been fully appreciated in the framework of phase-field models. 

It should be mentioned that a second obstacle for obtaining quantitative results on alloys 
is the strong contrast of the solute diffusivities between solid and liquid, which generates 
spurious solute trapping when the interface thickness is scaled up. This problem was solved 



by the so-cal 



approach 15, 



ed antitrapping current, which was introduced first in the coarse-graining 



21] , but has also been incorporated in the two-phase model |4j, |22] . Since the 



results of the present paper do not introduce major changes on this point, the results of Ref. 



15| will be taken over without a detailed discussion. After a change of variables, the model 



is almost identical to the model of Ref. 



15j . such that the asymptotic analysis developed 



there remains valid. This will be illustrated by numerical simulations that explicitly test the 



independence of simulation results of the interface thickness for the case of a lens-shaped 
phase diagram. 

The remainder of the paper is structured as follows. In Sec. [Til the standard phase-field 
model for the solidification of pure substances is reviewd for reference. Next, the grand 
potential formulation is introduced and motivated in Sec. IIII[ and illustrated by several 
examples in Sec. IIV1 a model with two parabolic free energy functions, a dilute alloy, and 
an alloy with a lens-shaped phase diagram. The relation of this model to other phase-field 
models is clarified in Sec. |V] and an example for numerical simulations is presented in Sec. 
IVIl Finally, the implications of the present findings for the further developments of alloy 
phase field models are discussed in Sec. IVIII 

II. SOLIDIFICATION OF A PURE SUBSTANCE 

The minimal model of solidification is considered, which implies some standard simplifi- 
cations: the densities of the solid and the liquid are taken to be equal, and heat transport 
is assumed to take place by diffusion only. As a consequence, no motion of matter needs to 
be considered, and the only transported extensive quantity is heat. 

Under these assumptions, the state of an inhomogeneous two-phase system can be com- 
pletely specified on a coarse-grained (mesoscopic) scale by two fields: a phase field <fi which 
indicates the local state (liquid or solid) of matter, and the internal energy density e. The 
conservation law corresponding to heat transport is the conservation of energy, which writes 



where k(4>) is the heat conductivity, which in general may depend on both variables, and T 
is the local temperature. 

If e is chosen as dynamic variable, the corresponding thermodynamic potential is the 
entropy. Consequently, the system is naturally described by an entropy functional 



d t e = -V • j e . 



(1) 



The heat current j e is given by Fourier's law, 




(2) 




(3) 



5 



where the local entropy density s depends on the fields as well as their gradients. Using the 
thermodynamic definition of the temperature, 

- = - (4) 
the evolution equations for the fields e and can be written in a variational form, 

d t <f> = (5) 

d t e = V ■ (k(0)Vt) = -V ■ ^(0)T 2 V^) . (6) 

The first equation expresses the local maximization of the entropy, which occurs at a rate 
given by the constant M^; the second is identical to the conservation law, Eq. ([I]). 

While this is a perfectly viable starting point which has been explored by several authors 
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24j . this formulation is rarely used in practice. Simulations are almost always carried 
out with models formulated in the variables <fi and T that can be obtained from free energy 
functionals. Besides historical reasons (the first phase-field models for solidification were 



formulated in this language 



25l . |26||). there are also formal considerations which make 
this approach preferrable. The main reason to choose the intensive variable T instead of 
the extensive variable e is that it directly controls the two-phase equilibrium, which makes 
it easier and more intuitive to identify the driving forces in the model. 

This point will be illustrated by obtaining the equations of motion from a free energy 
functional that is constructed using a purely phenomenological point of view. It will be 
shown below that the standard phase-field model of solidification is easily obtained as a 
special case. Let / S (T) and fi(T) be the free energy densities of pure solid and liquid, 
respectively, and let the corresponding equilibrium values of the phase field be <fi = ±1. The 
free energy is given by 

JF[<f>,T\ = [ /(</>, W,T)= / / int (0,V0)+^(0)/ s (T) + (l-^(0))/KT), (7) 
Jv Jv 

where the weighting function g s {4>) is given by 

9s{(p) = 2 ' ( 8 ) 

with g((f)) a function that satisfies g(±l) = 1 and g'(±l) = 0; hence, g s = 1 in the solid and 
g s = in the liquid. The term / int is given by 

/int = ^(V0) 2 + i// dw (0), (9) 



where a and H are constants of dimension energy per unit length and energy per unit 
volume, respectively, and /dw(0) is a double- well function with minima at = ±1. 

The motivations for this formulation are easily understood and are common to many 
phase-field models. The term f\ Qt creates domains where the phase field is close to its 
equilibrium values = ±1 (the minima of the double-well function), separated by diffuse 
interfaces. Therefore, far from the interfaces, the free energy density reduces to the one 
of the corresponding bulk phase. The term / int contributes to the free energy only inside 
the interfaces; this excess free energy represents the surface tension. The function g s ((f)) 
interpolates between the two free energy densities through the diffuse interface. 

A variation of the free energy functional with respect to the two fields and T yields 

Sjr = J {-^ 2 + #/dw(0) + ^ [/.CO - fi{T)\ } 50(f) 

+ [ufi^jp- + a - (io) 

where the prime stands for derivation with respect to 0. 

The variation of J 7 with respect to is the driving force for the phase transition. Two 
features are noteworthy: (i) the requirement that #'(±1) = ensures that the driving 
force vanishes outside of the interfacial regions, and (ii) since, at the melting temperature 
T = T m , f s (T m ) = fi(T m ), the "thermodynamical" part of the driving force is identically 
zero, independently of the value of 0. The latter property implies that the equilibrium 
interface profile in can be calculated from the term f int alone. This can be shown by 
seeking the equilibrium solution for a planar interface along the x direction, which can be 
obtained from the condition that the variation of T with repect to must vanish. At 
T = T m , this condition yields 

-ad xx <i> + Hf' dw {<l>) = (11) 

which implies that the solution of this equation is independent of the free energies f s (T) and 
fi(T). As a consequence, the surface free energy 7 (defined as the excess free energy due to 
the presence of the interface) is also independent of the bulk free energies and is given by 

7 = iV^H = IHW, (12) 

where J is a numerical constant that depends on the shape of the double well function / dw , 



7 



and the interface thickness W is defined by 

w = Vf (13) 

Therefore, the surface tension 7 and the interface thickness W can be freely chosen by 
appropriately fixing the two constants a and H, independently of the bulk properties. As 
stated in the introduction, the interface properties can thus be controlled independently of 
the bulk thermodynamics. 

The equation of motion for the phase field is obtained from the free energy functional by 
the standard variational procedure, 

d t <P = ~M^, (14) 

which expresses the fact that the system seeks to minimize its local free energy at a rate 
which is controlled by the constant M^. 

To obtain an evolution equation for the temperature field, the starting point is the ob- 
servation that, by definition, the variation of T with respect to T is equal to the negative 
of the local entropy density, 

s(T, 0) = = -g s {<t>)—^ (1 - 9s(W—^-- (15) 

Since / int was chosen independent of temperature, s is a local function of and T (no 
gradients are involved) and can be simply seen as the interpolation of the bulk entropy 
densities, 

s(<f>,T) = g s (<f>)s s (T) + (1 - gs^UiT). (16) 
The use of the thermodynamic identity de = Tds (valid at constant density) yields 

^m^r^^+^y). (it) 

Note that in writing down the second equality, it is assumed that there is no entropy produc- 
tion due to local dissipation, which is equivalent to the hypothesis that the transformations 
are reversible on the mesoscopic scale of a coarse-graining cell. Furthermore, the definition 
of the specific heat per unit volume is 

rur\ r ds ^ T ) TW , ^ 2 /,(T) n u , W) ( . 

C p (<f>,T) = T—^p — = -Tg s {(f))—^ T(l-g s (<p)) — . (18) 



The equations (IT7|) and (|T8|) can be combined with the energy conservation law, Eq. (JTJ), 
where the heat conductivity k now depends on the variables (f) and T, to yield 

ds -> r -» 

C p (0, T)^T = + V T) VTj , (19) 

which is the desired evolution equation for the temperature field. This equation can be 
further simplified by writing the heat conductivity as the product of the specific heat and 
the thermal diffusion coefficient Dt, and by using Eq. (fT6l) for the entropy density. The 
result is 



d t T 1 



C p (<f>,T)D T (<f>,T)VT +T[ Sl (T)-s s (T)}^-d t 4>}. (20) 



C P (0,T) 

Note that all quantities that appear in this equation except for the thermal diffusivity can 
be obtained from the bulk free energy densities f s {T) and fi(T). 

The free energy functional used in the standard formulation of the phase-field model 
can be obtained from Eq. ([7]) by linearizing the free energy density f(<p, V0, T) around the 
melting temperature T m , 



/(0,V0,T) = /(0,V0,T m )+^ 



(T — T m ) . (21) 



Using the definition of the latent heat per unit volume, L = T m [si(T m ) — s s (T m )}, as well 
as the fact that the free energies of solid and liquid are equal at T m , f s (T m ) = fi(T m ), this 
expansion yields 

T = / \a{V<Pf + Hf dw (<P) + ^r9(<P)(T - T m ), (22) 

JV * 

where the constant [s s (T m ) + si(T m )}/2 has been disregarded for simplicity. 

The equation of motion for the phase field is then obtained from this linearized functional 
by a variational derivative. To obtain an equation of motion for the temperature, it is 
usually assumed that the specific heat is independent of both temperature and the phase 
field, C p ((p,T) = C p . Then, it is easy to "guess" the correct equation by realizing that the 
latent heat released or consumed during the phase transformation appears as a source term 
in the diffusion equation for the temperature, 



d t T = V (D T (0,T)VT) + §-^-d t <f>. (23) 



This equation is indeed obtained from Eq. (120]) when a constant specific heat is inserted and 
the approximation T = T m is made in the second term on the right hand side. Note that, 
in contrast, the correct general form could not have been easily guessed from Eq. ( I23p . The 
underlying reason is that the linearized free energy functional formally yields a specific heat 
which is zero since all second derivatives with respect to the temperature vanish; therefore, 
thermodynamic consistency between the linearized functional, Eqs. ( 122]) . and the evolution 
equation for the temperature, Eq. (1231 . has been lost. The equations are nevertheless correct, 
since for the case of constant specific heat and constant latent heat, that is T[si(T) — s s (T) = 
T m [si(T m )— s s (T m )] = L, the internal energy density and the temperature are linearly related, 

e(0, T) = ^^t 6 ^ + C P (T - T m ) - 9 M L . (24) 



Then, Eq. ( 1231) can be directly obtained by combining the time derivative of Eq. (124]) with 
the energy conservation law, Eq. (CD). 



III. ISOTHERMAL ALLOY SOLIDIFICATION 

A binary alloy is a mixture of two pure substances A and B. For simplicity, it is assumed 
here that the atomic volume V a of both pure substances and of the mixture are all the 
same, and that hence the total number density of the alloy is a constant and equal in solid 
and liquid. Then, the only new field needed is the local composition (atomic fraction) c of 
"solute" (B) atoms. Furthermore, for constant atomic volume the chemical potentials of A 
and B atoms are not independent since removal of an A atom implies the addition of a B 
atom. This means that the only new intensive variable that needs to be considered is the 
chemical potential /i of the solute atoms. 



The starting point of the coarse-graining approach, as pioneered in Refs. [1CH12| , is a free 
energy functional that depends on the variables 0, T, and c. The chemical potential is then 
defined as the functional derivative of the free energy functional with respect to c. Since 
c is dimensionless, the chemical potential obtained by this procedure has the dimensions 
of energy per unit volume. This convention obscures the thermodynamic roles of the two 
variables: the relevant extensive variable from which a density should be defined is the 
number of B atoms. Therefore, the variable that is analogous to the internal energy density 
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in the pure substance model is the number density of B atoms, 

P = (25) 

where V a is the atomic volume (the constant volume occupied by one A or B atom). Then, 
the chemical potential defined by 

has the dimension of an energy, as is standard in basic thermodynamics. Since the nature 
of the variables is important for establishing the analogy between pure substance and alloy 
models, this convention for the chemical potential is adopted for the remainder of the paper. 
However, since it is customary to express free energy densities in terms of the composition 
rather than the number density, both p and c will be used in the following for ease of 
presentation, keeping in mind that the two variables are simply related by Eq. (|25|) . 
The number density is a conserved quantity, which implies 

d tP = -V-j p} (27) 

where j p is the mass current. For isothermal solidification, where the only thermodynamic 
driving force for mass diffusion is the gradient of the chemical potential, the mass current 
is given by 

j p = -M((f>,T,c)^fi, (28) 

where M(0, T, c) is the atomic mobility. Combining these equations and the definition of \i 
yields the equation of motion for p, 

dtp = V- (M(^T,c)V^ . (29) 

It will now be shown that this approach generally leads to a model in which bulk and 
interface properties do not decouple. To this end, it is useful to start again from the variation 
of the free energy functional, now in the variables and p. Since isothermal solidification 
is considered, there is no variation with respect to temperature. In order to simplify the 
notations, the variable T (which becomes a simple parameter for isothermal solidification) 
will be dropped from the free energy densities and the mobility from now on. The variation 
of T is 

&? = J j-^vV + HfM + *M [/s(c) _ /l(c)] \ 5m 

+ Va [qM^ + (1 - 9M)^] (30) 
11 



A crucial difference which the pure substance case is obvious: there is no simple argument 
which ensures that the "thermodynamic driving force" term proportional to / s (c) — fi(c) 
vanishes. Indeed, for two-phase equilibrium in an alloy, both the free energy density and the 
concentration vary across the interface. The value of all these quantities in the bulk phases 
at two-phase coexistence are obtained from two conditions: (i) the chemical potential must 
be the same in both phases, and (ii) the grand-potential density uj = f — pp must also be 
the same. Given the curves of free energy versus composition, the graphical interpretation 
of these two conditions is the well-known common tangent construction. 

Let us examine the consequence of these conditions for the phase-field model outlined 
above. Since the two variables - phase field and concentration - vary through the interface, 
the equilibrium interface profile is given by two coupled nonlinear differential equations. 
One is obtained from the condition of constant chemical potential, which remains valid in 
the diffuse interface picture, and reads 

^ = VagM?^ + V a (l - Q s m dJ ^ = ^(T) (31) 

where /i eq (T) is the equilibrium value obtained from the common tangent construction. This 
equation defines an implicit relation between the composition c and the phase field 0. 

The equation for the phase field, obtained as before from the condition that the variation 
of J 7 vanishes, is 

-ad xx <f> + Hf , dw (<f>) + ^- = 0, (32) 

where /th(0, c) = g s (<j>)fs(c) + (1 — 9s(4>))fi(c) denotes the "thermodynamic part" of the 
free energy density. Obviously, this equation becomes identical to Eq. (11 II) only if the third 
term is identically zero. The physical meaning of this condition can be made transparent by 
remarking that, since the concentration and the phase field are not independent variables 
any more under the constraint of Eq. f[3"Tj) . the variation of / t h with respect to <f), taking into 
account the constraint of Eq. ( 13T1) , is 



<?/th = d/th d/th dc = dfth p eq dc 

6</> ~ d<j) dc d(j)~ d<p V a d<p' [ ' 

Therefore, if 

% = 14 (Ah " ^ ^ °' (34) 

the equation for the equilibrium phase field profile reduces to Eq. (fTTT) ; in other words, the 

quantity f th — p eq p must be constant through the interface. Far from the interfaces, where 
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/i nt does not contribute, this quantity is equal to the grand potential density. Note that the 
common tangent construction implies that the two bulk values of the grand potential must 
be equal. However, for a general choice of free energy functions, there is no reason for this 
condition to be valid throughout the whole interface. As a consequence, the interface equation 
and all quantities that are obtained from its solution (surface tension, kinetic coefficients 
etc.) depend on the bulk free energy densities. For realistic values of the interface thickness, 
this dependence is small, but when the interface thickness is upscaled, large errors can occur. 

This fact has been recognized by several authors, and so far two different strategies have 
been followed to cure this problem. The first is to develop specifically designed free energy 
functionals that satisfy the condition of Eq. (1341) . but are valid only for certain choices of bulk 
free energies (see below). The second strategy is the one of the two-phase model [13|, |14| . 
in which two separate concentration fields, one for each phase, are used; the supplementary 
degree of freedom is then eliminated in such a way that Eq. AMI) is satisfied. 

The new idea put forward here is that a general solution to this problem can also be 
obtained in the coarse-graining spirit (using a single concentration field) when the model 
is derived from a grand-potential functional instead of a free energy functional. Indeed, 
the model formulated in the variables <p an d p is equivalent to the pure substance model 
formulated in terms of <fi and e: it has the same variational structure [compare Eqs. ([5]) and 
(|29p ], and both p and e are densities of extensive variables. To obtain the equivalent of the 
more successful pure substance models formulated in the variables (f) and T, alloy models 
should be formulated in the variables and p; the corresponding thermodynamic potential 
is the grand potential. 

grand-potential functionals have been used extensively in classical density functional the- 
ory (see [23] for a review). In the context of phase-field models, a grand potential functional 



has 1m -on introduced to study solidification with density change 28|, I29]. However, in all 
the cited works, the density is retained as the fundamental field that is used to evaluate 
the functional. In contrast, the grand potential in its role as a thermodynamic potential 
depends on the chemical potential. If the goal is to have a complete formulation of the 
problem in terms of the dynamical variable p, the grand-potential functional Q should be a 
functional of the field p. In thermodynamic equilibrium, this field is just a constant which 
is equal to the thermodynamic equilibrium chemical potential, but out of equilibrium p can 
depend on space and time. Therefore, the field p that appears in the free energy density 
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needs to be eliminated in favor of p. This is simple if the values of p and p are related by 
a local and invertible function. A free energy functional of the form of Eq. (jTJ), taken with 
free energy densities that depend on T and c, is a good starting point since it contains no 
nonlocal terms in p such as (Vp) 2 . Moreover, for functions f s (c) and //(c) that are convex 
in c, the relation between p and c is monotonous and hence invertible. Thus, it is possible 
to switch from c to p as the dynamic field. After this operation, the number density is not 
a fundamental free field any more, but is obtained as a local functional derivative of the 
grand-potential functional with respect to the local chemical potential, 

p =yjMiA, (35 ) 

op 

Note that the above requirements (no square gradient terms in c and convex free energy 
functions) implies that the present method cannot be applied to systems that exhibit phase 
separation. 

In analogy with Eq. (JTJ), the grand-potential functional is 



tt[(j),p] = I L)((j), V0,/i) 

Wmt(0, V0) + g s ((f))u s (p) + (1 - g s ((j)))ui(p), (36) 



v 

where Ui nt is identical to fi nt , and the grand potential densities of the bulk phases are 
obtained by a Legendre transform of the free energies, 

u„(p) = f»{c) - pp (u = s, I). (37) 

This procedure can be easily performed for any convex free energy function, either analyti- 
cally or numerically. Note that this transformation implicitly uses the equivalence between 
statistical ensembles (canonical and grand canonical) on the mesoscopic scale. This is con- 
sistent with the general philosophy of the coarse-graining approach, which assumes that 
thermodynamic quantities can be defined on the scale of a coarse-graining cell. 
The variation of the grand-potential functional (at constant temperature) is 

sn = J |-<rvV + HfU4>) + [".(/<) - wiM]} 

+ (38) 

Since, at solid-liquid coexistence, u s = u>i, the equilibrium interface equation obtained from 
the condition of vanishing variation with respect to 4> is identical to Eq. ffTTj) . as desired. 
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Furthermore, using Eq. (|35|) . the variation of Q with respect to p yields 
p{(j),p) ■■ 



5Q dus(fi) , dui^j) 



(39) 



which is the equivalent of Eq. f fl5|) for the entropy density. Note that this can also be 
rewritten as 

p(0,T,/i) = gMps(T,fi) + (1 - g s {4>)) Pl {T,p) (40) 
with p v = du u /dp [y = s,l). It is useful to restate this equation in terms of c for future use, 



c(0, p) = V a p((f), fi) = g s ((p)c s (fi) + (1 - g s {4>))ci{fi), 



(41) 



where obviously c u (p) = V a du u /dp. Two-phase coexistence is characterized by a constant 
chemical potential, p = p cq (T); the corresponding composition profile through a solid-liquid 
interface is 

c eq (0) = c(0,/i cq (T)) = g s (<j>)C(T) + (1 - ^(0))cf(T). (42) 

The equations of motion for and fi are now formulated following the same steps as for 
the pure substance model. The phase field evolves toward a minimum of the grand potential, 

5Q 



-M A 



Ma 



(43) 



this equation shows that the thermodynamic driving force for the phase transition is the 
difference in grand potential densities. The evolution equation for the chemical potential is 
obtained by taking the time derivative of Eq. fj4*0|) . which yields 



/ dp((j),fi) dp((f),p) 



(44) 



It is useful to define the quantity 



X\<PiW = — a = 9sK<f>)—R — + (1-^(0))- 



(45) 



op op op 

which will play a role similar to the specific heat in the pure substance model. The symbol % 
is chosen here because this quantity can be seen as a generalized susceptibility 30j . Further- 
more, the mobility can be written as the product of x(0, p) and a solute diffusion coefficient 
D((p,p). Indeed, in the bulk phases, for monotonous (and hence invertible) functions p Sy i(p), 

dpM 1 1 



dp 



dp v {p)ldp V a 2 dy v {c)/dc^ 



(46) 
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which is the well-known thermodynamic factor (Darken factor 



The combination of 



these definitions and Eq. (j44|) with the mass conservation law, Eq. (j27|) yields 

dt ^ = W^{^' [ D ^»)X(<!>,»)V»] -^p-{ps{»)-pi(»))d t( f>y (47) 



an equation completely equivalent to Eq. (1201) for the temperature in the pure substance 
model. 

IV. EXAMPLES 

The model is completely specified for any set of free energy functions for solid and liquid 
by the definitions of the grand-potential functional, Eqs. ( )36|) and (1371) . and the evolution 
equations ( )43|) and ( ]47l) . In order to illustrate some of its properties, it is useful to work 
out several explicit examples. First, it is shown that the equivalent of the linearized pure 
substance model is obtained from parabolic free energies. Next, it will be shown that the 



dilute alloy model of Ref. 



15] can be recovered using this formalism. Finally, the more 



general case of an ideal solution model will be treated. 
A. Parabolic free energies 

The simplest phenomenological approximation for free energy functions for fixed equilib- 
rium compositions c e s q and cf 1 at some temperature T are two parabolas, 

Uc)= l -e u {c-C)\ {v = s,l) (48) 

where e s and e/ are constants with dimension energy per unit volume. The chemical potential 
in each phase is 

M=^ = Wc-4 (49) 
which can of course be inverted to yield c as a function of fi in each phase, 

+C (50) 



V a e 

The grand potential densities are then obtained from the Legendre transform, uj v = /„ — fip, 
where Eq. (|50|) is used to switch variables from c to ji, 

= -\t£- - f C- (51) 
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Of course, the use of the definition p = —du/dp together with c = V a p yields again Eq. (I50|) . 
The equilibrium chemical potential for two-phase coexistence is obtained by the condition 
cu s (fi eq ) = Ui(fi eq ). The solution ^ cq = corresponds to the common tangent between the 
bottoms of the parabolic wells, and the equilibrium compositions are equal to c e s q and cf 1 . 
For e s ^ ei, a second solution exists which corresponds to a common tangent that is tilted 
and yields different values for the equilibrium compositions; this solution is not of interest 
here. Since p cq = 0, we have /i = p — fi eq , that is, in this model p is directly the deviation 
from the equilibrium value of the chemical potential. 

Inserting these grand potential densities in Eqs. fl4Tl) and (T4"51) yields the expressions for 
c and x as a function of and p, 



c(4>,fi) 



1 -g.(4>) + ^(1-^(0)) 



A* + c«0.(0+cf(l-0.(0) (52) 



X(<f>, H) = 771^(0) + 772-0- - 9s{<t>)). (53) 



V a e s V a 6 t 

It is easy to see that when e s = q = e, the resulting model is equivalent to the standard 
pure substance model. Indeed, in this case the difference between the compositions in the 
two phases is independent of /i and identical to Ac = c^ q — c e s q , and the susceptibility is just 
a constant, x = ^-/(V a 2 e). This corresponds to the approximations of constant latent heat 
and constant specific heat, respectively. It is easily verified that Eq. (l5"2l) becomes identical 
to Eq. (I2"lj) . with Ac and x rep 



in several phase-field models 



acing L and C p , respectively. This analogy has been used 
18N20|. The present formalism makes it possible to generalize 



this model and to use e s ^ €[, that is, parabolas with different curvatures. 



B. Dilute alloy 



The phase-field model for a dilute alloy of Ref. 15j starts from the free energy densities 



U(T, c) = tf(T) + e„c + %^ (cln c - c) , [y = s, I) (54) 

where ff(T) and fi A (T) are the free energy densities of pure A, e s and q are again constants 
with dimension energy per unit volume, ks is Boltzmann's constant, and the last term in 
Eq. (|54"|) is the dilute limit of the entropy of mixing. 

The calculation of the grand potential densities is straightforward and yields 

eA k B T ( p - V a e h 



^_^exp (55) 
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The compositions as a function of the chemical potential are given by 

ffi-V a e v \ . . 

C » = eXP {-J^)- (56) 

It is obvious that these compositions satisfy, for any value of the chemical potential, the 
partition relation 

c s = kci (57) 

with the partition coefficient k given by 

k = ^( Va{ ?V a) )- (58) 



k B T 

The equilibrium chemical potential for a given temperature is again obtained from the 
condition u s (T,fi cq ) = u>i(T, /i eq ) which yields 



A/rr\ { A/rr\ "'B 



/.CO - fftT) 



k B T 



V n 



a 



exp, ^^r exp v k R r 



(59) 



The left hand side can be expanded for temperatures close to the melting temperature T m ; 
the right hand side can be rewritten in terms of q q = q(/i cq ) and the partition coefficient. 
The result is 

^ L (r-T ra ) = ^ f f(fc-1). (60) 

If, on the right hand side of this equation, T is approximated by T m , which is justified in 
the limit c < 1, the standard dilute alloy phase diagram is obtained, 

T — T 

cf = (61) 

where m = k B T^(l — k)/(V a L) is the liquidus slope. 

When these expressions are now used in the complete grand-potential functional, the 
composition as a function of <fi and /x becomes 

c(0, = gM) exp (^J^) + (1 " 9s(<f>)) exp (^^) • (62) 

A specific feature of the dilute alloy model is that this expression can be factorized into 
two functions that depend only on <p and /i, respectively. Moreover, the latter can again be 
rewritten in terms of the liquid composition, which yields 



C (0,/i) = Q(/x) [1 - (1 - % s (0)] . 
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(63) 



As a result, the expression for the susceptibility is also quite simple and reads 

1 ( U — Vatl \ j \ n C(0, fl) 



' 1^1 1 (1 k)g t 



(64) 



V a k B T ^\k B Tj l y V a k B T 

Inserting these expressions in Eqs. fj4~3]) and ff4T|) leads to model equations identical to 
those of Ref . 15] . The difference of the grand potential densities is 



Va 



H - V a e s \ ( fj, - V a e t 

exp — - — — — — exp 



(65) 



k B T J \ k B T 

The free energy difference can be expressed in terms of the equilibrium chemical potential 
using Eq. (1591) . and after some algebra one obtains 

k B T 



Va 



-cf (1 - k) 



exp 



knT 



- 1 



(66) 



As in Ref. [15j, two dimensionless variables are now introduced. The first, 

A* A'-cq 



k B T 



(67) 



is the dimensionless deviation of the chemical potential from its equilibrium value. This 
variable can also be expressed in terms of the composition and the equilibrium composition 
at two-phase coexistence using the definition of Eq. ( 1421) as 



The second dimensionless variable, 



u = In 



U 



Ceq(0) ' 

e u - 1 



(68) 



(69) 



is a dimensionless supersaturation. This can be seen by inserting the identity e u 
c((fi, fi)/c eq ((j)) obtained from the preceding equation, which yields 

c((j),ii) - C eq (0) 



U 



(70) 



(1-A;)c cq (0) • 

When the driving force (the grand potential difference) is expressed in either of these vari 



15| is obtained. 



ables, the same evolution equation for the phase field as in Ref. 

Matters are slightly more complicated for the second evolution equation. Ref. 15] uses 
an evolution equation for the composition (or, equivalently, for the dimensionless super- 
saturation U) rather than for the chemical potential. For the dilute alloy model, such an 
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equation can be obtained starting from the general evolution equation of the chemical po- 
tential, Eq. ( |4"T|) . or (more simply) from the original formulation of the mass conservation 
law, Eq. (|2"7|) . For this purpose, the chemical potential has to be eliminated in favor of c 
or U. This is possible because Eqs. (16"7"|) and Eq. ( 1ST?]) can be inverted and combined with 
Eq. (JM]) to yield 

H = /icq + k B T In -^j- (71) 

c / t 1 " I 1 - k )9s{<i>)\ 

and 

c = c^l - (1 -%,(</>)] [1 + (1-W (72) 

A straightforward calculation then yields the variational form (without the antitrapping 
current) of the evolution equation for the composition of Ref. [15|. These steps will be 
discussed in more details later on; here, it is important to stress that the possibility to 
switch from an evolution equation for fi to one for c by an exact transformation is specific 
to the dilute alloy model: this procedure only works because the function c(0, fj) is easily 
inverted, which is not the case in more general models, as will be seen below. 



C. Ideal solution model 



In an ideal solution model, the free energy is a weighted average of the pure substance 
free energies of A and B plus an entropy of mixing term, 

f„(T, c) = (1 - C )/*(T) + cf*(T) + ^ [cine + (1 - c) ln(l - c)] . (73) 

* a 



Contact with the notations of the previous examples can be made by setting t u = f„ — 
the free energy becomes 

/„(T, c) = tf{T) + ce v {T) + ^ [cln c + (1 - c) ln(l - c)] . (74) 

* a 



The chemical potential is 



dfu c 

fj, = — - = V a e v + k B T In . (75) 

op 1 — c 



This can be inverted to yield the conentration in each phase as a function of /a, 

c v (n) = 7 — - — r-. (76) 



1 + exp (^) 
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The grand potential densities are 



= ft + 



ln(l - c) = 




In 1 + exp 



) 



(77) 



As before, these expressions can now be used to define the interpolated composition and the 
susceptibility, 



where the latter has been expressed in terms of the functions c Sj /(/x) because this leads to a 
simpler expression. 

The equations of motion for the ideal solution model are obtained by inserting these 
expressions in the general evolution equations, Eqs. ( )43|) and ( )47l) . The equilibrium chemical 
potential and the phase diagram can be calculated analytically, but the resulting expressions 
are quite complicated and not of interest here. The important point here is that for this 
model it is not possible to transform the evolution equation for \i into one for the composition 
c: whereas the functions c„(/x) for the composition in each phase as a function of [i can be 
easily inverted, the same is not true for the interpolated composition given by Eq. (1781 . As 
a consequence, it is easy to obtain c from \x for given and T, but hard to obtain [i from c. 

V. RELATIONS WITH OTHER PHASE-FIELD MODELS 

A. Equivalence to the Kim-Kim-Suzuki model 

I L 

In the two-phase model [14| . the interface region is seen as a phenomenological super- 
position of the two bulk phases, with a weighting function h s ((f>) that interpolates between 
liquid (h s = 0) and solid (h s = 1). The main difference to the coarse-graining approach is 
that the two-phase model starts with two separate composition fields for the solid and the 
liquid, c s and q. The free energy density and the "true" composition (in the coarse-graining 
sense) are then written as 




(78) 



X(</>,A*) 



c s (/i)(l - c s (fj,))g a (<f)) + c; (//)(! - ci(fj))(l - g s ((p)) 

V a k B T 



(79) 



f(cf>, c) = h s (4>)f s (c s ) + (l - h s {4>))fi{ Cl ), 

21 



(80) 



c = h a (<f>)c a + (l-h,{<f>))c l . (81) 

Since there are two equations but three variables (0, c s , and q), an additional condition is 
needed to close the system: the chemical potentials of the two coexisting phases are required 
to be the same, 

= df a (<t>,c a ) = dfi((j),ci) 
dc s dci 

Taking this implicit relationship between c s and q into account, the number of independent 
fields is reduced to two, and two evolution equations for, say, and c can be written down. 

This approach is completely equivalent to the grand-potential formalism outlined above, 
with the difference that the dynamical variable is c instead of \x. To see this, set h s {4>) = 
g s (<t>)- Equation (15T1) becomes identical to Eq. (|41|) . Furthermore, since the compositions 
c s and q in Eq. (HT1) are defined as functions of the variable /i, Eq. (1521) is automatically 
satisfied. Finally, since all the equations of Ref. [ijj] are developed by analogy with the pure 
substance model, it is not surprising that the evolution equations, Eqs. (31) and (32) of 
Ref. [u| are identical to the evolution equations (H31 and (H7|) of the present paper once all 
the notations have been properly translated: the driving force for the phase transformation 
is the difference in grand potential density and the quantity f cc that appears in the evolution 
equation for the concentration in Ref. 14J is identical to l/(V a 2 x) here. 

In summary, the two-phase model can be obtained in a fully variational manner from 
a grand-potential functional. Note that, in this point of view, the introduction of two 
separate composition fields is not necessary any more: the fundamental dynamical field is 
the chemical potential, and the two fields c s and q are simply obtained as the derivatives of 
the bulk grand potential densities with respect to fi, whereas the "real" composition at any 
space point is given by the functional derivative of the grand-potential functional. 



B. Relation to the phenomenological two-phase model 

At this point, it is useful to discuss the respective merits of the grand-potential and the 
two-phase approaches. At first glance, the former seems to be more advantageous: the num- 
ber of fundamental fields is equal to two (as compared to three for the two-phase model), and 
whereas fi has to be obtained from c s and q in the two-phase model by solving Eq. (1521 . c s 
and q are obtained from \x by simple derivatives in the grand-potential formalism. However, 
an analysis in terms of computation time reveals that matters are not so simple. For the 
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sake of concreteness, consider the simplest case of constant but distinct solute diffusivities 
in the solid and the liquid. In this case, the evolution equation for the composition reduces, 
in the bulk phases, to the simple linear diffusion equation. In contrast, the evolution equa- 
tion for jj, has several nonlinearities (even in the bulk) due to the presence of the factors 
x{4>i /-0 > which are in general nonlinear functions of p. This means that the numerical effort 
to integrate the equation for \i in the bulk is higher that the one for integrating an equation 
for c. In the interfaces, the grand potential formulation does not require much additional 
effort, whereas the nontrivial Eq. (1821) has to be solved in the two-phase model. However, 
since the interface regions usually represent only a small fraction of all the grid points in a 
numerical simulation, the equation of motion for p may not always be advantageous from a 
computational point of view. 

The computational disadvantage of the grand-potential formalism in the bulk can be 
alleviated by an additional step, which also brings to light the direct relation of this approach 
to the phenomenological two-phase model developed by the Access group [13]. The idea is to 
make a change of variables and to replace the chemical potential field by another continous 
field that plays the same role. In order to obtain the standard diffusion equation for this 
field in the liquid phase, the appropriate field is the density (or the composition) in the 
liquid, pi or q. Indeed, under the hypothesis that p = dfi/dp is an invertible function of 
p, the function pi = dui/dp is just its inverse function, according to the properties of the 
Legendre transform. Furthermore, p s = du s /dp is also a function of p and can therefore 
expressed as a function of pi, 

p s {pi) = ps{p{pi))- (83) 
Eliminating p in favor of pi in Eq. fH7|) yields 

d -- & {* ■ h< «> + ™ ^ - * 3 A ■ (84 » 

where xi — d 2 ui/dp 2 = is the susceptibility of the liquid phase. Unsing the fact 

that x{4>i p) = Xs(/-O<?s(0) + ~ 9s(4>)) by definition, this equation can be rewritten as 

1 



d t pi 



Xi(pi)[l - (1 - g s ((j))Xs/xi)] 



g'{4>) 



+ 



D{^ Pl )xi{pi)[l - {l-9s{4>)XslxiWpi 
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For a dilute alloy, p s /pi = k and Xs/xi — k according to Eqs. ( 158|) and (|64j) . and the above 
model becomes identical to the one of Ref. 13| . The grand-potential formalism thus allows 
for a generalization of this model to arbitrary phase diagrams. 



C. Local supersaturation approximation 



Further progress can be made by a simple approximation. The crucial point is the 
relation between composition and chemical potential: its nonlinearity penalizes the grand- 
potential formalism in the bulk and makes the resolution of Eq. ( 182]) in the two-phase model 
non-trivial. The idea is thus to replace the exact relation between chemical potential and 
composition in the interfaces by an approximate one that will make it possible to write a 
simple equation for the composition. This approximation will be called local supersaturation 
approximation: it exploits the fact that, for slow solidification, the chemical potential in 
the interfaces is close to the value for two-phase equilibrium. This suggests to use a Taylor 
expansion of the composition around the equilibrium composition profile in the interface, 



dc 

c-c cq (0)=- 



- A*eq) = VaX(4>, ^eq)0" _ A*eq) 



(86) 



which can be easily inverted to yield 



Ceq(0) 



(87) 



Furthermore, an expansion of the driving force around the equilibrium chemical potential 
yields 



U a (fJL) - Ui(jj) 



„eq _ eq 



V, 



-(fj.-fj.eq). 



(88) 



Inserting these expression in the evolution equation for the phase field and the mass conser- 
vation law yields a simple set of equations for and c, 



1 



Or 



TTf, C-C cq (0) 

^idw 7T- Ac 



Va X(0,/^eq) 



d t c = V 



D(0) X («^ eq )V (^ eq(0) 



• /^eqj 



Note that the second equation can also be rewritten after applying the chain rule as 



d t c = VL>(0)Vc + V 



D(4>) 



Ac + 



Xi 



eq 



X(0,A*eq) / 



V0 



(89) 
(90) 

(91) 
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where xTi = xi^i A'eq)- This latter form displays explicitly the two driving forces for solute 
diffusion that establish the equilibrium solute profile and that lead to solute redistribu- 
tion out of equilibrium: composition gradients and differences between the thermodynamic 
potentials in the two phases. 

Before proceeding further, it is useful to relate the quantities that appear in the above 
equations to the phase diagram of the binary alloy, characterized by the curves c| q (T) and 
q q (T), or equivalently by the coexistence line fi eq (T). The quantities xl q are related to the 
liquidus and solidus slopes, m v = dT/dd^J 1 , 



1 dd? _ dc 
m v dT dfi 



VaX?^- (92) 



Moreover, the quantity dfi eq /dT can be evaluated using a Clausius-Clapeyron relation for 
the fi-T coexistence line, 

d/ieq _ L _ LV a 



Combining these results yields 



dT TAp TAc 

XZ = ( 94 ) 
V a Lm Sy i 

Thus, the susceptibilities are inversely proportional to the liquidus and solidus slopes, and 
therefore x e s q /Xi q = m i/ m s- Defining an effective partition coefficient by the ratio of the 
liquidus solpes k m = mi/m s , the evolution equations of the phase- field model can be further 
simplified. Indeed, the susceptibility along the equilibrium profile is 

x(0,/i eq ) = X c s q 9s(<P) + X/ q (l - 9si4>)) = XHI " (1 " U^(0)], (95) 

and the evolution equations for and c become 

1 a a V72, ufi 9'{4>) Ac c-c cq (0) 

^ =lV '" Hf " ~ — « i - (i - *-)«.(« m 



d t c 



D{(j>)[l-{l-k m )g t 



c - c eq (0) 



(97) 



l - (l - km)g s {4>). 

These equations are very similar to the ones of the dilute alloy model, except that the parti- 
tion coefficient k has been replaced by the effective partition coefficient k m which depends on 
the temperature. In that sense, this approach bears some similarity with the method used 
in Ref. [32] to construct a quantitative phase-field model for arbitrary phase diagrams: the 
free energy curves are first approximated by a dilute alloy phase diagram with "effective" 
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(temperature-dependent) partition coefficient, melting temperature, and liquidus slope; the 
equations of motion for the dilute alloy model are then applied with these effective param- 
eters. Here, the approximation is directly in the evolution equations, and can be applied 
in a straightforward manner for arbitrary free energy functions. Also note that the two ap- 
promixations are not completely equivalent. For instance, the effective partition coefficient 
in Ref. [32| is defined by the ratio of the compositions, c^ q /q q , which is equal to the ratio 



mi / m s only for a dilute alloy. 

It should also be noted that this approximation is not equivalent to the approximation 
of constant concentration gap Ac (equivalent to parallel liquidus and solidus slopes) used 



in Refs. 



18 
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33). This can be easily seen from Eq. (j86p by considering a constant 



chemical potential deviation 5fi = // — /i eq (generated, for example, by a local curvature of 
the interface): the shifts in concentration on the two sides of the interface are proportional 
to xTii respectively, and hence inversely proportional to the liquidus slopes, as they should 
be. 

The above approximation has been called "local" for two reasons. First, this emphasizes 
the fact that the approximation of the relationship between composition and chemical poten- 
tial is needed only in the interfaces, while it leaves the bulk evolution equations unchanged. 
Second, it is anticipated that this method should be applicable to situations in which the 
temperature field varies slowly with time and over large length scales, such as in directional 
solidification or in thermosolutal models with realistic values of the Lewis number. For this 
purpose, it should be sufficient to apply the above equations with the local value of the tem- 
perature for each point of the interface, which implies that Ac, xTi anc ^ are n °t constants 
but vary between different interface points. 



D. Non- variational model and antitrapping current 

The next step is to incorporate into the model two features that have been widely used 
to increase the precision and performance of phase-field models. The first is motivated 
by the fact that a non-variational model can be more efficient for computational purposes 
than a strictly variational formulation. This was first highlighted by Karma and Rappel jlj], 
and their method has since been used in many other models. In the terms of the present 
formulation, it amounts to keeping the interpolation function g(<j)) in the evolution equation 
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of the phase field, but using a different interpolation function for the number density and 
the susceptibility. Let h(<p) be a function that has the property h(±l) = ±1, and let 

*.<#) = (98) 
Then, the concentration is defined as 

p(0,aO = = pMKW + pi(aO(i - M0) (99) 

» a 

instead of Eq. fj40l ; as a consequence, the equilibrium composition profile given by Eq. fj42|) 
is also modified and becomes 

c eq (0) = c^h s (<j>) + cf (1 - h s {(f>)). (100) 

The susceptibility is still defined as the derivative of the number density with respect to the 
chemical potential and becomes 

X (0,T,/i) = gp( y /i) = h s (<j>) Xs (T,p) + (1 - h s (<j>)) X i(T,p). (101) 

The second feature is the so-called antitrapping current, which was developed in order 



to counterbalance spurious solute trapping [15|, |2l| . Indeed, if the solute diffusivity in the 



solid is substantially lower than in the liquid, as is usually the case in alloy solidification, the 
upscaling of the interface thickness magnifies the solute trapping effect, whose magnitude 



is proportional to the interface thickness 34). To restore local equilibrium at the interface, 



as appropriate for low-speed solidification, an additional solute current is introduced which 
"pushes" the solute out of the freezing material, and which is given by 

Jat = aW[pi(fM) - p s (fx)]nd t (j), (102) 

where n = — V0/|V0| is the unit normal vector pointing from the solid to the liquid, W is 
the interface thickness, and a > is a constant to be determined by a matched asymptotic 
analysis [lsj. The current thus defined is proportional to the interface velocity (via the 
factor d t (p) and to the composition difference between solid and liquid, and is directed from 
the solid to the liquid for a solidifying interface, for which d t <p > 0. The mass conservation 
law, Eq. (!27j) . reads now 

dtp = -V (J P + Jat) = V (M(0, p)Vp - aW[ Pl {p) - p s {p)}nd t <p) ■ (103) 
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Taking these modifications into account, the evolution equation for the chemical potential, 
Eq. (JIT]) , is replaced by 



dtH 



1 



D{(j), A*) V/i - aW[pi(/i) - p s (p)]hd t <f) 



(104) 



the evolution equation for remains unchanged. 

In the local supersaturation approximation, the composition difference in the expression 



for the antitrapping current is approximated by its equilibrium value, pi — p s 



c e s q )/V a = Ac/V a . The evolution equation for the concentration, Eq. f[9"T]) . is then replaced 
by 



dtc= V 



D(0)[l-(l-A^)/i,(0)]V 



c - C eq (0) 



1 - (1 - fc m )fc s (0) 



(105) 



E. Relation to the quantitative dilute alloy model 

As usual, the parameters of the phase-field model have to be related to the quantities 
that appear in the sharp-interface theories by matched asymptotic analysis. The complete 
asymptotic analysis for the general Eqs. ( )4"3]) and ( l4Tj) will be presented elsewhere. Here, 
only the behavior of the model in the local supersaturation approximation will be discussed, 



because for its analysis the analogy to the dilute alloy model 15[ can be exploited. 

In order to apply directly the results of Ref . [15| , it is useful to cast the evolution equations 
in dimensionless form. From Eqs. ( 189]) and (190]) . it is clear that the dimensionless variable 
that generalizes the quantity U of the dilute alloy model defined by Eq. ( 1691) is 



U 



cq 



Ceq(</>) 



Ceq(0) 



X(0,/ieq) AC 



Ac[l 



(106) 



Indeed, this expression can be rewritten using the fact that c eq (4>) = q q [l — (1 — k)h. 
where k = c^ q /q q is the standard partition coefficient, 



U 



k)h a (<f>) c-c eq (0) 



1 - k 1 - (1 - k m )h s (4>) c eq (</>) 
which reduces to Eq. (170]) for a dilute alloy, since k = k m in this case. 



(107) 
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In terms of this variable, Eq. ( 190]) becomes 

[1 - (1 - k m )h.{<f>)]d t U = V{£>(0)[1 - (1 - k m )h s {(j))]VU 

+ aWh[l + (1 - k m )U]d t (j)} 
+ [! + (!- k m )U\d t h s {4>). (108) 



A form formally identical to Eq. (69) in Ref. [15| is obtained by choosing a particular 
interpolation for the diffusion coefficient D((p), namely by setting 

D{<f>)[l - (1 - k m )h s {(j>)] = D iq {<j>), (109) 

where Di is the solute diffusivity in the liquid, supposed to be constant. Note that, since 
the left hand side is actually the product of the diffusivity and the susceptibility, strictly 
speaking the function g(0) is an interpolation of the mobility rather than of the diffusion 
coefficient. The final result for the evolution equation is 

[1 - (1 - k m )h s {4>)\d t U = V (p iq {4>)VU + aWn[l + (1 - k m )U\d t ^j 

+ [l + (l-k m )U]d t hM- (no) 

This equation is indeed identical to the one used in the asymptotics of the dilute alloy model, 
except that the dilute alloy partition coefficient k is replaced by the ratio of the solidus and 
liquidus slopes k m . This is very natural, since this quantity controls how the composition 
difference between solid and liquid depends on the chemical potential at the interface. 
The replacement of c by U in the evolution equation for the phase field leads to 

_U = ,v*0-^-^P>. (in) 

2 VaXi 

This equation is now divided by the constant H, which amounts to non-dimensionalizing 
the free energy and grand potential densities (since H has dimension of energy per unit 
volume) . Furthermore, the function g is now chosen to be the standard fifth-order polynomial 
g((p) = 15(0 — 20 3 /3 + 5 /5)/8, the function h((p) = <fi, and the double-well function to be 
/dw = 1/4 — 2 /2 + 4 /4. The resulting equation reads 



r<9 t = jy 2 V 2 + 0-0 3 - (1-0 2 ) 2 A£7, (112) 



where r = 1/ (M^H) is the relaxation time for the phase field, W = \fo~j~H. is the interface 
thickness defined by Eq. ffl3|) . and the constant A is given by 

15 (Ac) 2 

16 HV a \T 



A = ^T7=fe- (H3) 
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Equation (I112p is identical to the standard evolution equation for the phase field 

As announced previously, the results of the asymptotic analysis of Ref. jjjj] can now be 
exploited since Eqs. (IllOp and (11121) are identically to the model analyzed in this reference. 
Therefore, the variable U obeys, in the liquid, the free boundary problem 

d t U = AV 2 f/, (114) 

= -doJC - fiV n , (115) 
[1 + (1 - k m )U int ]V n = -D^UU, (116) 

where K, and V n are the local curvature and interface velocity, respectively, do is the capillary 
length, and (3 the kinetic coefficient. Equation f 1 1 1 5 j) is the generalized Gibbs-Thomson 
equation, and Eq. (11 16)) is the Stefan boundary condition that describes mass conservation 
at the phase boudary. 

In terms of the phase-field parameters, the capillary length and the kinetic coefficient are 
given by 

W 

do = a x — (117) 



T 



with ai = 5\/2/8 and 02 = 0.6267; these values are identical to those obtained by Karma 
and Rappel 

In terms of physical quantites, this expression for the capillary length is in fact identical 



to the standard thermodynamic definition 30|. Indeed, the number a\ = 5\/2/8 quoted 



above is equal to (15/16)/, where /, the numerical constant defined in Eq. (1121) . is equal to 
2\/2/3 for the standard fourth-order double well potential used here. With the help of these 
relations, Eq. (11171) can be rewritten as 

IWHV^y^ 



(120) 



Then, the use of Eqs. (1121 and (J46j) yields 



do 



(Ac) 2 (Ac)2 &m 
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VI. NUMERICAL TESTS 



The relations found in the preceding section are now used to perform some illustrative 
simulations on a concrete model system. The ideal solution model is chosen, with the same 
parameters as those used in Ref. 12[ to model the Nickel-Copper alloy. This alloy exhibits 
a typical lens-shape phase diagram with a rather narrow coexistence zone. Concretely, the 
free energy densities given by Eq. f J73|) are used; the free energy differences (between solid 
and liquid) of the pure substances are given by 

jNi 

f m {T) _ f m {T) (T _ t»% (121) 

m 
tCu 

f Cu {T) _ f cu {T) = (T _ T m^ (122) 

m 

with the melting temperatures = 1728 K and T^ u = 1358 K and the latent heats 
L Nl = 2350 J/cm 3 and L Cu = 1728 J/cm 3 . The molar volume is taken to be 7.42 cm 3 , and 
the surface tension 7 = 3.3 x 10~ 5 J/cm 2 . 

For a temperature of 1710 K, the equilibrium concentrations are c e s q = 0.045988 and 
q q = 0.058098, the partition coefficient is k = 0.7916, the ratio of the liquidus slopes 
is k m = 0.8017, and the capillary length calculated by Eq. (fl20|) is d = 6,426 x 10~ 6 
cm. Isothermal dendritic solidification is simulated in two dimensions. The anisotropy 
needed to obtain stable dendritic growth is introduced in the standard way [l[ 21 j by 



making the surface tension dependent on the angle 9 between the interface normal and a 
crystallographic axis, here chosen to coincide with the x direction. A standard fourfold 
anisotropy, 7(0) = 7(1 + e4Cos(46*)) is used, with 64 = 0.025. The initial composition of the 
liquid is chosen as (c e s q + c e l q )/2, which corresponds to a supersaturation of 0.5. 

Four simulations are carried out with the equations of motion in the local supersatura- 
tion approximation, with A = 1.596,3.192,4.788, and 6.384, which corresponds to interface 
widths of W = 116,232,348, and 464 nm. For each simulation, the relaxation time of the 
phase-field equation is chosen such as to eliminate interface kinetics (/3 = for all orien- 
tations). The steady-state growth velocity of the dendrites is measured as described in 
Ref. and the result is displayed in Fig. [TJ It exhibits the behavior that is typical for 
quantitative phase-field models: the simulation results converge to a constant value with 
decreasing interface thickness W, and the convergence is roughly quadratic in W, which can 
be expected since all terms linear in W have been eliminated by the model formulation. It 
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FIG. 1. Growth velocity of a two-dimensional dendrite in the Ni-Cu alloy at 1710 K, for different 
choices of the interface width W. 

is not surprising to find such behavior, since the model used here is essentially identical to 



15 



2l|. 



the dilute alloy model for which rapid convergence with W has been demonstrated 
Moreover, the parameters chosen here are in a region of the phase diagram in which the 
dilute approximation should still work quite well. However, this is not a limitation of the 
approach: simulations at T = 1410 K, where the ratio of the liquidus slopes k m = 1.2609 
is very different from the partition coefficient k = 0.9590, yield a similar convergence plot. 
This shows that the model can be applied to alloys with arbitrary phase diagram. 



VII. SUMMARY AND PERSPECTIVES 



The most important conclusions of the present work can be summarized as follows. 

1. A phase-field model for alloy solidification has been obtained in a completely varia- 
tional framework, starting from a phenomenological grand-potential functional that 
is a simple sum of bulk and interface contributions. In this model, the two dynamic 
variables are the phase field and the chemical potential field. A complete analogy 
can be established with the standard phase-field model for the solidification of a pure 
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substance, in which the variables are the phase field and the temperature field. The 
main difference between the two sets of equations is that the relation between the 
composition and the chemical potential (the extensive and the intensive variable) is 
nonlinear, whereas the relation between temperature and internal energy is usually 
assumed to be linear. 



2. The resulting model is shown to be completely equivalent to the model of Kim,Kim 



and Suzuki 



14j , with the chemical potential replacing the composition as the dynamic 



field. As a result, in the present model no "partitioning" of the solute into coexisting 
phases is necessary; the (auxiliary) compositions in each of the phases can be simply 
obtained from the chemical potential. With an additional change of variables, the 

n 

phenomenological Access model [13J can also be recovered and extended to general 
alloy phase diagrams. These developments show that these two-phase models, despite 
a seemingly quite different starting point, can in fact also be obtained by a coarse- 
graining procedure if the appropriate thermodynamic potential is used. 



3. The equations of the new model can also be written in terms of a phase field and the 
composition field. They thus have essentially the same computational complexity as 



the original phase-field models for alloy solidification |10L [11] derived in the coarse- 
graining framework. However, in contrast to the latter models, they have a decisive 
property which is required for quantitative simulations: bulk and interface thermody- 
namic properties can be adjusted independently. This difference in behavior is due to 
the different interpolations of the relevant thermodynamic potentials (free energies in 



Refs. 



10 



11] , grand potentials here). 



4. With an additional approximation - a linearization of the relation between chemical 
potential and composition inside the interfaces - the model becomes equivalent to 



the quantitative dilute alloy model studied in Refs. [15l . |2l|. This feature makes it 
possible to include the additional antitrapping current as in that model, and to apply 
the detailed asymptotic analysis of Ref . |15[ • As demonstrated for one particular case 
here, efficient and accurate simulations are thus possible for arbitrary alloy phase 
diagrams. 
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Numerous interesting perpectives for future work arise from the present results. First, 
the model has been worked out here for isothermal solidification only, but it can be easily 
extended to non-isothermal situations: coupled equations for the phase field, the tempera- 
ture field and the chemical potential field can be developed by following the same steps as 
done here for each transport field separately, taking also into account off-diagonal elements 
in the Onsager matrix of transport coefficients as well as cross-derivatives of the thermody- 
namicpotentials. This is anticipated to yield a generalization of the thermosolutal model of 



Ref. 



16|. 



Moreover, the generalization of the formalism to mult i- component systems should be 
straightforward. This is particularly interesting because general models for multi-component 



multi-phase solidification have been developed in the two-phase framework [4J, |22|, |35| . In 
these models, the determination of the compositions in the individual phases for given global 
composition requires to solve a generalization of Eq. f )82|) . which represents a set of coupled 
nonlinear equations (one for each component). A formulation in terms of the chemical 
potential would completely avoid this problem and thus potentially offer important gains in 
computational performance. 

It should also be mentioned that the two-phase approach has been used in other contexts. 



One example is the treatment of fluid flow in solidification 36[ and two-phase flows 371. |38| . 



where two separate velocity fields, one for each phase, are introduced, in contrast to "coarse 



graining" models in which a single velocity field is used 39[. It would be interesting to 



reassess the relations between these different models in the light of the present findings. 
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